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Abstract 

An effective control strategy with a non-homogeneous soil profile for horizontal Ground Source Heat Pumps (GSHPs) 
was investigated in this paper. Steps toward development of a comprehensive model to consider the effects of an 
intermediate blanket overlaying the ground pipes were described. The model incorporates the effects of a variety 
of surface energy fluxes to provide an accurate estimate of the ground thermal regime. The developed model was 
utilized successfully in conjunction with the Genetic Algorithm (GA) evolutionary search algorithm to obtain the 
optimized operational parameters for a GSHP in three different climate conditions. A properly sized and engineered 
non-homogeneous soil profile demonstrated the potential to boost the capacity of GSHP systems to a significant level. 
The potential benefits of a recycled product, Tire Derived Aggregate (TDA), as an insulation blanket was assessed via 
the optimization algorithm. TDA was demonstrated to be more effective in the heating mode in a cold climate (Buffalo) 
by increasing the energy extraction rates from the ground by about approximately 15 % annually. TDA's effectiveness 
was less pronounced in a relatively moderate climate (Dallas) for heating purposes with efficiency improvements of 
about 4 % annually. The annual percentage increase in efficiency for cooling season with TDA blanket was higher for 
Buffalo (7.6 %) compared to Dallas (3.5 %). For the cooling only purposes (Miami), a high conductive intermediate 
layer (saturated clay) exhibited greater potential to enhance the efficiency in coldest months in a warm climate. The 
results are highly suggestive of the beneficial application of a layered system to increase the performance of GSHPs. 
A shift in design approach toward consideration of more control strategies for the ground pipe side of GSHPs is 
suggested based on the model results. A field demonstration project with TDA blanket has been constructed at the 
University at Buffalo to further investigate the findings with regard to field conditions and other characteristics of TDA 
material. 

Keywords: Ground source heat pump, tire derived aggregate, non-homogeneous, optimization, genetic algorithm, 
energy efficiency, control 



1. Introduction 

Every ground source heat pump (GSHP) consists of two main parts; the ground side and the heat pump inside 
the building. The ground side, also referred to as the heat source/sink, consists of the ground pipe network and the 



"Corresponding author's address: 204 Jarvis Hall, University at Buffalo, Buffalo, NY, 14260. Tel: (716) 380-5255. Fax: (716) 645-3667 
Email addresses: ar92@buf f alo . edu (Amir Rezaei-Bazkiaei), ehsandehObuf f alo . edu (Ehsan Dehghan-Niri), 
gdargushObuf f alo . edu (Gary F. Dargush), sweberSbuf f alo . edu (A. Scott Weber) 

'Department of Civil , Structural & Environmental Engineering, University at Buffalo, Buffalo, NY-14260 
department of Mechanical & Aerospace Engineering, University at Buffalo, Buffalo, NY-14260 

3 Smart Structures Research Lab, Department of Civil, Structural & Environmental Engineering, University at Buffalo, NY 14260 



Preprint submitted to Energy and Buildings 



Monday 9 th April, 2012 



circulating pump. The heat pump unit itself consists of sub-units which handle the thermodynamic relationships 
between the working fluid and the load side (indoor air or indoor water network); compressor, condenser, evaporator, 
the expansion device and/or blower fans or circulating indoor water pumps. The working fluid flows through the 
ground pipes and exchanges heat with the surrounding soil medium where it gains/loses heat in heating/cooling 
modes respectively. The fluid at the outlet of the pipe enters the heat pump where the thermodynamic cycle of 
the heat exchange between the working fluid and the refrigerant is responsible for heat delivery or dissipation in 
heating and cooling modes respectively. A successful design of GSHP involves careful selection and sizing of both 
the ground side and the indoor unit so that heating/cooling loads are met year-round with least waste heat/electricity 
inventory. This selection process usually involves employment of simplifying engineering assumptions. For example, 
one of the most well-known ground pipe sizing semi-empirical formulas is based on the concept of thermal resistance 
calculation for the soil, pipe and working fluid [1], Resistance factor calculation is based on the line source theory 
which assumes a constant source of heat propagates continuously through the soil medium of interest. An assumption 
which is not a real representative of the ground pipe operation but is assumed to be within reasonable range of error 
for an engineering practice. 

The heat pump sizing starts with the careful estimation of the heating/cooling loads for the prospect building via 
accepted methods such as the bin method. The bin method is based on a full survey of the total number of hours in 
a design year that a certain air temperature range (bin) occurs. The heating/cooling budget analysis forms the base 
for the peak load GSHP design. It is noteworthy that methods like bin method do not take into account the sequence 
or real-time occurrence of the peak load demand in time. In other words, the question is the duration of a worst 
case load condition rather than when and for how long the load sustain. Therefore, it does not matter in such design 
procedures what is the climate and/or source condition before and after the worst case condition, a matter that can 
influence the overall performance of the system. Moreover, these methods usually assume a constant source/sink 
strength (ground temperature) based on the maximum and minimum observed data or semi-experimental formulas 
[1] which again does not reflect on the dynamic variation of the boundary conditions and the time dependent nature 
of the system's functionality. Inherent tendency of the traditional design methods like this to stray from the path of 
optimized operation has made the dynamic, real-time analysis of GSHPs the focus of many recent research works 
[2, 3]. The following design steps after the heat budget analysis involve ground pipe sizing and heat pump selection. 

One of the most essential parts of a GSHP design, after a meticulous selection of the pipe length and size, is the 
proper selection of the working fluid temperature range that enters the ground. International ground source heat pump 
association (IGSHPA) suggests an equation mainly based on experiment (equation 1) which relates the inlet water 
temperature to the heat pump to maximum, minimum and average ambient air temperatures. It is the responsibility 
of the designer to select the appropriate optimal entering water temperature to the heat pump that guarantees efficient 
performance. Given the complexities of an accurate building load demand evaluation (users consumption habits, level 
of building insulation, changes in occupancy, and specially sudden climatic changes which are far different from the 
long term statistical design year data, etc.) and its causal relationship with the GSHP system units, the optimization 
of the overall performance of a GSHP becomes a real challenge. Add to the problem, the ground characteristics 
variation in different regions to make it even more challenging to track the dynamic building load demand and keep 
the heating/conditioning process as optimized as possible. This paper aims to help designers of GSHPs achieve 
optimal design specifications by utilization of an evolutionary optimization algorithm. The central motive behind 
this research is to obtain optimized values for the functional parameters of a GSHP system which can be applied 
practically via employment of control measures such as variable refrigerant flow heat pumps. 
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With the ever increasing need for preserving energy, the idea of utilization of capacity control strategies for GSHPs 
was initiated by researchers about two decades ago [4]. The capacity control practices in their core deal with the most 
feasible engineering options to dynamically meet the building load requirements and minimize the waste energy in 
the heating/conditioning process. As Madani summarizes [5], capacity control usually involves either the control on 
the components of a GSHP system (e.g. the compressor, condenser, ground pipes, etc.) or a change in design con- 
figurations for different seasons or advanced control algorithms. The complex nature of the interplay of relationships 
between the units of a GSHP is in such a way that changes in operational parameters of one unit can affect the overall 
performance of the system. From the modeling point of view, researchers have tried to tackle these complexities by 
dynamically modeling the interactions between these units via packaged softwares such as TRNSYS [6, 7, 8, 9, 10, 2] 
or in combination with thermodynamic data base models such as EES [3] or by self developed programs [11] which 
are capable of modeling the dynamic between these units. The general structure of such modeling efforts often con- 
sists of separate first-level models (e.g. ground pipe and heat pump) which each have sub-models of second-level (e.g. 
compressor, evaporator, etc.) with higher details. The functional relationship between the sub-models is constructed 
via operational parameters such as water/brine and refrigerant flow rates. The level of complexity of these sub-models 
significantly depends on the purpose of the modeling effort as Madani [5] describes. 

Despite numerous efforts to impose control measures on the heat pump side of a GSHP system, there has been 
put little thoughts into the potential of having control on the ground characteristics. This is perhaps because of a 
predisposition to believe that ground related works usually are associated with extra capital investment that makes 
any modification less favorable than making changes to the heat pump unit. Our findings [12] suggest that utilization 
of a layered system above the pipe burial depth might be worth more attention. The layered system demonstrated 
the potential to increase the energy extraction rates from the ground to about 17 %, in the peak heating month, 
compared to the homogeneous case with the local soil profile. This paper aims at assessing the effectiveness of a 
non-homogeneous soil profile for different climate conditions by obtaining the optimized functional parameters for 
these different climates. 

Our focus in this paper is mainly on the analysis of the system efficiency and optimization via control on the 
source/sink (ground) side with an eye on the heat pump unit. A simple semi-empirical equation (1) for the entering 
water temperature to the pump was adopted from the international ground source heat pump association (IGSHPA) 
design manual [1] to make the link between the ground pipe and the heat pump. The ground pipe model, on the other 
hand, is an elaborate version of a surface energy balance model which unlike thermal resistance method is capable of 
solving for the temperature distribution of the entire soil profile in order to obtain the outlet water temperature from 
the ground pipes. A detailed description of the surface energy balance equations, parameters and the solution methods 
is presented in [12]. A summary of the heat fluxes is listed in the 'surface energies' section. Our motivation behind 
development of such ground model arises from the desire to investigate the potential benefits of a non-homogeneous 
soil profile on the ground pipe performance. The developed model was employed to investigate the effectiveness of 
utilization of a recycled product of tire industry as an insulation material above the burial depth of the ground pipes. 
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where: EWTh'. entering water temperature in heating mode, EWT C : entering water temperature in cooling mode, 
EWT min : minimum design entering water temperature, EWT mean : average design entering water temperature, EWT max : 
maximum design entering water temperature, T airmean : average annual ambient air temperature, T airmin : minimum an- 
nual ambient air temperature, T airmax : maximum annual ambient air temperature, T air : air temperature at the time of 
simulation 

Tire Derived Aggregate (TDA) is the proposed insulation material to be used as the intermediate layer in the non- 
homogeneous model. TDA mainly consists of chopped pieces of used tires in a variety of nominal sizes ranging from 
1 to above 5 inches [13]. The idea of production of TDA from used tires, also referred to as tire chips, tire shreds, 
and tire mulch, was first initiated by Humphrey [14]. New applications for TDA, mostly in civil engineering, were 
proposed based on its unique physical characteristics. TDA's relatively low density compared to conventional backfill 
makes it a viable alternative fill material where a lighter fill material is desired in construction works [15]. Its low 
thermal conductivity on the other hand, provoked the thought of utilizing TDA as an alternative insulation material 
in road base insulation and some agricultural applications to modulate temperature fluctuation on the soil surface 
[14, 13, 16]. There are only a few studies in the literature that focus on the properties of TDA measuring physical 
properties [16], compaction densities for different sizes of chips [13] and also thermal properties of tire chip samples 
[13, 17, 18]. A list of thermal properties and density of tire shreds reported in these studies is presented in Table 1. 
A thermal conductivity value of 0.29 Wm~ lo C~ 1 , specific heat of 500 JKg lo C 1 and density of 720 KgmT 3 were 
chosen as representative properties of TDA material for the purpose of the analyses in this paper. 



Table 1 : TDA thermal and physical properties in literature 



Reference 


Material 


Thermal conductivity 


Specific heat 


Density 






(Wm- lo C- 1 ) 




(Kgm- 3 ) 


Wappet (2006) 


Tire shred 


0.564 


507 


641 


Shao(1995) 


Tire chips 


0.149-0.164 


NA 


513 


Humphrey (2002) 


Tire chips 


0.29 


1.15 


720 


Moo-young (2003) 


Tire chips 


NA 


NA 


1060-1100 



In the following sections, the development of a model for horizontal ground pipes in a non-homogeneous soil 
profile was described first. The model was then used in conjunction with the developed genetic algorithm optimization 
model to search for the optimized operational parameters of the GSHP. The optimization was carried out for same 
pipe characteristics and different cities across United States to evaluate the effects of climatic conditions on the select 
optimized parameters. Detailed description of the methods and results are presented in the following sections. 

2. Material and methods 

A brief summary of the surface energy balance equations are first introduced. The second subsection describes the 
ground pipe specifications and also the way pipes are modeled from the physical domain into the numerical domain. 
Numerical solution of the ground pipe model and the parameter settings for the genetic algorithm are described in the 
following sections. 

2.1. Surface energies 

The surface boundary condition takes into account the effects of energy balance due to variety of mechanisms 
responsible for surface-ambient heat interaction. The total energy balance on the ground surface (Q,, Watts) can be 
written as [19, 20]: 
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Qt = Qc + Qe + Qh + Qie + Qu + Qsi + Qp 



(2) 



Q c - Conduction heat flux through snow layer or ground surface (Watts) 

Q e - Turbulent exchange of latent heat (Watts) 

Qh- Turbulent exchange of sensible heat (Watts) 

Qi e = Emitted long wave radiation heat flux (Watts) 

Qu= Incoming long wave radiation (Watts) 

Q si - Solar radiation reaching the surface of earth (Watts) 

Q p - Heat flux due to precipitation (Watts) 

For a detailed description of each heat flux input parameters and formulas refer to [12]. 
2.2. Ground pipe modeling 

The ground pipes configuration consists of a horizontal pipe of total length of 61 meters (L) and | inch diameter, 
3 m distance between inlet and outlet, buried at the depth of 2 m (j P i pe ), with the working fluid flow rate of 3 gallons 
per minute (0.19 kg/s). Assuming there is no thermal interaction between pipes, the solution domain (Figure 1) was 
considered to be from the center line of the pipe to the mid-span of the distance between pipes (1.5 m), in the x- 
direction, and from ground surface to the farfield (5 m) in the y-direction. Soil thermal conductivity (K s ) and thermal 
diffusivity (a s ) were selected based on the soil and rock classification guideline [21] to be equal to 1.67 Wm~ l °K~ l 
and 66 • 10 -8 m 2 s~ x , respectively. For the sake of comparison, this physical setting of the ground pipe was assumed to 
be the same for other cities in this analysis. 

To account for the three-dimensional behavior of the pipe and the surrounding soil, the effect of the working fluid 
flow rate was considered along the pipe direction. The third dimension of the problem was modeled by splitting the 
physical domain in the pipe direction into a series of cross sections (slices) of the soil profile for each time step, 
including the nodal temperature of the fluid at the pipe's location. Figure 2 shows how the slices are spaced in the pipe 
direction to cover the temperature distribution of the entire 3D domain. At each time step, the nodal temperatures of 
each cross section were obtained and subsequently updated for the next slice along the pipe's length via equation (4) 
to achieve a temperature distribution of soil and fluid at the end of the pipe (L). The same process was repeated for 
the next time steps until the end of the simulation time. A schematic representation of the cross sections is depicted 
in Figure 1. Detailed description of the numerical schemes and results are described in the following sections. 
3. Numerical algorithm 

The three-dimensional temperature distribution in the soil media was modeled by solving the governing heat 
conduction equation and incorporating the heat flow rate forced by the circulating water. The temperature gradient 
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Figure 1: Configuration of pipe, TDA and soil layers, and solution domain discretization [12] 




Figure 2: Schematic of slices of 3D domain in pipe directions [12] 



in the pipe material is small enough to be neglected allowing the heat equation to be solved in a two-dimensional 
geometry by inclusion of the fluid temperature in the domain [22]. The governing equation (3) was solved for the 
entire domain in two-dimensions for the pure heat conduction first, and the solution was then updated along the pipe's 
length to obtain the solution in three dimensions. 



ldT(x,y,t) d 2 T(x,y,t) d 2 T(x,y,t) 



a dt 



dx 2 



dy 2 



(3) 
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The output water temperature equation (4) along the pipe direction, calculated based on the analytical solution for 
the energy balance between surrounding soil medium and pipe [23], constructs the link between the fluid and the soil 
temperature in the model. 

T f ,ou, = T s - (T, - T fJ )e*£ (4) 



where a is the thermal diffusivity of the medium through which heat travels, Tf t0Ut is the fluid temperature exiting a 
pipe of length L (m). T s is the surrounding soil temperature, Tfj is the initial water temperature entering the pipe, and 
K s is the soil thermal conductivity (Wm"'°C _1 ). m is the mass flow rate (Kgs^ 1 ) and C p j is the specific heat of the 
working fluid (Jkg~ l °C~ l ). Boundary and initial conditions for the solution domain are described as follows. Initial 
temperature distribution of the soil profile (r,) was obtained from the Kusuda model, equation (5)[24]: 

Ti = T(x,y),t = 0 
— fl x — x 

fl x — u, a — Amax 

f x =0,x = 0 

Q,(W/m 2 ),y = 0 

T(y, t) = T(y, t) Kusuda , y = y max 



T(y, t) Kusuda = T avg + T amp e y cos(2tt-^ -y (5) 

where T avg is the average annual surface temperature, T amp is the amplitude of fluctuation of the annual surface tem- 
perature, y is the depth from ground surface (m), a s is soil thermal diffusivity (m 2 s~ l ), P is the duration of a year in 
seconds, and t is the time of the year in seconds. 

Partial derivation of temperature in x-direction was written with a central difference scheme: 
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where: 



$ 2 X T = Ti-ij - 2T U + T i+Uj 
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(7) 
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Figure 3: Two dimensional grid of solution domain [12] 



Using the same discretization for y-direction, the governing heat equation therefore can be expressed in terms of nodal 
temperature values in its explicit form as in equation (8). Next, a scheme was selected to march through time. As 
depicted in Figure 3, the solution domain was discretized in n x + 1 nodes in x and n y + 1 nodes in the y-direction, of 
which inner domain was used to solve for the temperature distribution of each cross section of the physical domain. 
The nodes on the boundary were separated to force boundary conditions in x and y direction. 
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A fully explicit finite difference solution scheme was employed to solve for temperature distribution of the so- 
lution domain. The model was constructed in MATLAB. Time steps of 1800 seconds and space discretization of 
Ax = Ay = 0.1 m were chosen based on a stability analysis of the model undertaken in the author's previous work 
[12]. 

Fully explicit finite difference formulation of the heat conduction equation takes the form of equation (9). 

'lu 1 r(Tl hj - 27?. + + K^y-i - 27£ + Tf J+l ) (9) 

where: 

aAt aAt 

r x = r y = r (10) 



x ~ (Ax) 2 ' y ~ (Ay) 2 ' x ~ y 
Rearranged form of the explicit equation will take the form of equation (1 1): 
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r +1 = rTU- + rTf +1> . + rT^ + r^ +1 + (1 - 4r)7£ (11) 

Once the homogeneous case formulation was performed, the next step was adjustments to the difference equations 
to take into account the effects of the internal boundary conditions on the top and bottom of the intermediate layer. 
To meet the temperature and flux conditions on the internal boundaries, the energy balance on the interfacial nodes 
of intermediate layer and soil was obtained by writing the nodal fluxes and the energy storage term for the control 
volume around each node. Figure 4 depicts the energy balance for a node on the top and bottom of the intermediate 
layer interface, the control volume around the node, and the heat fluxes involved. The resulting energy balance is 
shown in the equation (12): 

dTij 

qi + #3 - qz ~ #4 = (pC p ) avg ■ Ax ■ &y—^-, (12) 

where: 
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Figure 4: The intermediate (TDA) layer top and bottom interfaces energy balance [12] 

and the values of heat fluxes are: 



T n — T n 

91 = y 



(13) 
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Substituting the flux equations into the energy balance equation and rearranging the equations results in the introduc- 
tion of new parameters in the difference equation as described below. 



Kav g At K TDA At K s At 

r> r T DA = zr^?> r s = ~r-^\ TlTTo ( 17 ) 



avs (pC p ) avg (Ax) 2 ' lUA (pC p ) avg (Ax) 2 ' s { P C p ) avg (Ax) 2 

The difference formulas obtained previously for the explicit method were revised to account for the heat flux 
through the TDA top and bottom layers resulting in the equations (18) and (19). 
TDA top layer: 

T n+1 = r mg Tl h . + r avg T? +l . + r s + r T T? j+l + (1 - 2r avg - r s - r T )T n i4 (18) 



TDA bottom layer: 

T" +1 = r avg Tl uj + r avg Tf +l j + r T T n ihX + r s Tf j+l + (1 - 2r avg - r s - r^T^ (19) 
4. Genetic algorithm 

The idea of utilization of Genetic Algorithm (GA) in engineering applications was inspired by natural selection. 
This selection represents the concept of the higher survival chance of the fittest individual in its environment through 
successive generations to find the optimal design parameters among the others. More information about GA optimiza- 
tion can be found in [25]. The main advantages of GA over traditional optimization algorithms are that it does not 
require and does not depend on gradient information of the objective function; instead it uses a population of design 
points and randomly utilizes information from each generation to the subsequent one; and, there is a potential to make 
the population converge to the Pareto-optimal set. The latter property is a crucial appeal of GA in multi-objective 
optimization problems that enables us to use GA in combination with a Pareto-set filter [26] to obtain a near approxi- 
mation of the entire set of non-dominated solutions. The specifications of the GA implementation in the present study 
are listed in Table 2. 



Table 2: Gaoptimset options used in the analysis 



PopulationSize 


500 


CreationFcn 


@ gacreationuniform 


Generations 


50 


FitnessScalingFcn 


@fitscalingprop 


PopulationType 


doubleVector 


SelectionFcn 


@ selectiontournament 


EliteCount 


2 


CrossoverFcn 


@crossoversinglepoint 


CrossoverFraction 


0.7 


MutationFcn 


( @ mutationuniform, .02) 


MigrationDirection 


forward 


Display 


'iter' 


Migrationlnterval 


15 


UseParallel 


'always' 


MigrationFraction 


0.1 







Given the complex nature of the relationships between a GSHP unit, designers often should heavily rely on their 
level of experience to select the best heat pump design. Even with a huge deal of experience, these designs are 
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potentially distant from the required heating loads of a building due to the common practice peak load design criteria. 
Several researchers in the heat transfer literature have come to the conclusion that optimization schemes can benefit 
the design process where the effects of the climatic conditions is to be closely considered in an efficient design process 
[27, 28, 29, 30, 31, 32]. In their comprehensive work on vertical [30] and horizontal [31] ground source heat pumps, 
Sanaye modeled the thermodynamic cycle of the heat pump in conjunction with the thermal resistance pipe model to 
obtain the optimized operational parameters. Sayyaadi does a similar study via exergy analysis of a vertical ground 
source heat pump [32]. Authors of this paper are specifically concerned with the effects of the climatic conditions on 
the selection of the optimized operational parameters for a horizontal GSHP in the presence of a non-homogeneous 
soil profile. The optimization algorithm is supposed to provide the designers of the system with beneficial information 
regarding the ground pipe design where the absence of the experience can cause a wide deviation from the optimized 
cost and operation path. 

To evaluate the potential benefits of the intermediate layer in different climate conditions, three cities across the 
nation were selected in order to perform the optimization analysis. Buffalo's climate condition requires space heating 
for majority of the year (550 cooling-degree-days), so Buffalo was assumed the heating dominated city with eight 
months heating (starting on October) to be representative of the cold climate. Miami case (4383 cooling-degree- 
days) was analyzed for pure conditioning to represent the warm climate. Simulation for Dallas was done with the 
assumption of six month heating (November-April with cooling-degree-days less than 150) to be the representative 
of a mild climate condition. To simplify the introduction of the weather data, the air and dew-pint temperature, and 
solar radiation values for these cities were introduced to the model via estimation of these inputs by cosine functions 
(equations 20, 21, and 22). Weather data were obtained form the National Operational Hydrologic Remote Sensing 
Center (NOHRSC) and curve-fitting was undertaken in Microsoft Excel environment to obtain the input parameters 
to the cosine models. A summary of the input parameters to the model for these cities is presented in table 3. 



Table 3: Weather data for selected cities as input to the model 



City 


Buffalo 


Dallas 


Miami 


Ta-avg ( C) 


9.3 


19.7 


23.7 


T a-amp ( C) 


14.1 


12.2 


6.8 


a a i r (month) 


10.1 


10.1 


10.3 


Tdew-avg ( C) 


4.3 


9.4 


18.1 


T dew-amp ( C) 


12.2 


11.6 


7.3 


®dew {month) 


10.2 


10 


10.4 


Q si . avg (Wm- 2 ) 


142 


168 


182 


Qsi-amp (Wm' 2 ) 


75 


59 


33 


a Qsj (month) 


9.2 


10.2 


10.1 


U s (ms~ l ) 


4.2 


4.1 


4 



t oc ■ 

Tair = Ta-avg + T a-amp COS(2tT- - -jy) (20) 

where T a -avg is the average annual air temperature, T a -amp is the amplitude of fluctuation of annual air temperature, 
and a air is the fitted cosine model phase difference for air temperature, calculated based on the start time of the 
modeling on October first (the assumed heating season start time). 
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T dew point — 1 dew-avg + T dew-amp COS(27T^ 12 ^ 

where T dew-avg is the average annual dew-point temperature, T dew-amp is the amplitude of annual dew-point tempera- 
ture variation, and ctdew is the fitted cosine model phase difference for dew-point temperature. 

t Uq, 

Qsi = Qsi-avg + Qsi-amp COS(27T- - — ) (22) 

where Q S i- avg is the average annual solar radiation reaching the surface, Q S i-amp is the amplitude of fluctuation of the 
solar radiation throughout the year, and the a Qs . is the fitted cosine model phase difference for radiation. 

A genetic algorithm (GA) optimization scheme was designed to obtain the optimum values of the intermediate 
layer configuration (thickness and position), working fluid type, and inlet fluid temperature ranges. The results of 
the simulation will provide insights on the benefits of introduction of an intermediate layer on the performance of a 
GSHP throughout the year. The model outputs will clarify what are the optimal properties and configuration of the 
intermediate layer and what are the optimal entering water temperatures to the ground throughout the year to achieve 
the maximal efficiency. The ultimate motivation behind the analysis is to determine whether the operating parameters 
for different months of the year are in a range that confirms the benefits of utilization of a capacity control strategy or 
a set of selected operational parameters can be used for all months without deviating from the optimized condition. 



Table 4: Working fluid properties 


Index 


Fluid 


Thermal 


Specific heat 


Density 


Dynamic 


number 




conductivity 






viscosity 






(Wm- lo C- 1 ) 


(JKg- l °C- 1 ) 


(Kgm- 3 ) 


(Kgm- l s- 1 )- 10~ 3 


1 


Water 


0.6 


4183 


998.3 


1 


2 


6 % Propylene-glycol and water 


0.476 


4140 


1010 


1.5 


3 


13 % Propylene-glycol and water 


0.432 


4100 


1010 


1.9 


4 


18 % Propylene-glycol and water 


0.408 


4060 


1020 


3 


5 


24 % Propylene-glycol and water 


0.389 


4020 


1020 


6.3 



Seven input variables (decision variables) were chosen to feed the core finite difference model in each run of 
the GA. These inputs comprise working fluid properties, minimum, mean and maximum entering water temperature 
values, the intermediate layer thickness, position and thermal properties. To let the evolutionary algorithm search in 
a broader spectrum of potential answers, the algorithm was allowed to choose between a range of practical working 
fluid properties (table 4, extracted from the IGSHPA guideline [1]) and also a range of common soil properties in 
ground heat pump works (table 5, extracted from the manual for the soil and rock classification for design of ground 
heat pumps [21]). 

The values of heat extraction/dissipation rates were calculated based on equation (23) for each run after selection 
of these parameters. The time averaged values of outlet (T" l u f) and inlet (T™ g ) working fluid temperatures for the 
simulation period (monthly or seasonal) are used to calculate energy extraction/dissipation rates in heating/cooling 
modes. The heat pump work was calculated based on the friction factor and Reynolds number described in [33]. The 
objective function used in the simulation is equal to the reciprocal of difference between the energy extraction rate and 
the circulating pump energy consumption rate. A single-objective genetic algorithm model was used subsequently to 
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search for the parameters that minimize the objective function or, in other words, maximize the net energy extrac- 
tion/dissipation rates. It was assumed that maximum heating/cooling energy demand from the ground pipes does not 
exceed 1.5 Kw in any of the cities while the optimization was aimed towards maximizing the extraction/dissipation 
rates. A population size of 500 and 50 generation was selected for each run of the model. A total number of 32 pro- 
cessors of 48 GB memory were used for each monthly run of the model which resulted in runtime of about 7 hours. 
A total number of 64 processors of 48 GB memory were used for each annual run which resulted in the runtime of 
about 65 hours for each year. 

Energy = m-C p , r (TZf-T™ 8 ) (23) 



Table 5: Intermediate layer thermal properties 


Index 


Intermediate 


Thermal 


Thermal 


number 


material 


conductivity 


diffusivity 






(Wra- lo r') ( 


cm V)10- 8 


1 


TDA 


0.29 


58 


2 


Sand 


0.77 


45 


3 


Clay 


1.11 


54 


4 


Loam 


0.91 


49 


5 


Saturated silt or clay 


1.67 


66 


6 


Saturated sand 


2.5 


93 



3 3 




Days (starting on October 1st) 




0 50 100 150 200 250 300 350 



Days (starting on October 1st) 



Figure 5 : Annual air and solar radiation variation in the selected cities 



5. Results and discussion 

The optimized monthly parameters for three cities are presented in tables 6 to 8. The main motive behind doing the 
optimization separately for each month was to have an estimate of what the range of the optimal operating parameters 
are and how the optimal parameters vary from month to month. The selection of the intermediate layer material type 
by the algorithm assures a non-homogeneous profile that provides the best performance. It seems from the results 
that TDA blanket demonstrates more benefits in certain climate conditions and within each climate also its potential 
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advantage is more highlighted in some months more than others. Results for Buffalo (table 6) show that TDA is 
selected as the dominant intermediate layer for eight months of the year of which five months are in heating season. 
In the cooling season, TDA was not selected only in July where the saturated clay tends to exhibit a dominant effect 
as an intermediate layer. 



Table 6: GA results for Buffalo 



Month 


Fluid 


E W T m i n E W T max 

\ <- > \ L ) 


EWT mean 


Thickness 


Position 


Intermediate 

llltlLClltll 


Energy 
f Watt "I 

^ VV O.LL } 


% 


Oct 


3 


4.0 


5.2 


0.4 


1 A 

1.4 


Saturated clay 


ceo 

558 


6.5 


Nov 


1 


-3.0 


5.2 


0.7 


0.6 


TDA 


492 


4.7 


Dec 


2 


-3.0 


5.5 


0.9 


0.5 


TDA 


515 


13.2 


Jan 


2 


-3.0 


8.0 


1.0 


0.5 


TDA 


474 


18.8 


Feb 


2 


-3.0 


6.8 


1.0 


0.7 


TDA 


375 


19.4 


Mar 


1 


-3.0 


5.2 


0.5 


1.0 


TDA 


252 


6.8 


Apr 


3 


-2.6 


5.2 


0.8 


1.0 


Saturated clay 


179 


9.1 


May 


3 


4.7 


5.2 


0.4 


1.3 


Saturated Clay 


295 


11.7 


Jun 


1 


40.4 


35.8 


0.5 


0.5 


TDA 


1500 


2.5 


Jul 


2 


42.4 


32.5 


0.3 


1.4 


Saturated Clay 


1500 


1.2 


Aug 


4 


43.0 


41.6 


0.3 


1.0 


TDA 


1500 


1.2 


Sep 


4 


45.6 


34.5 


0.4 


0.5 


TDA 


1500 


1.3 


Avg. Heating 
Avg. Cooling 




-1 

42.8 


5.9 
36.1 


0.7 
0.4 


0.8 
0.8 




393 
1500 





The results for Dallas (table 7) on the other hand show a different trend where TDA has been selected as the 
dominant intermediate layer mostly in warmer months of the year whereas the colder months where saturated clay 
yields the optimal efficiency. The reason behind this observation probably lies in the difference in the phase change 
angle of the annual temperature and solar radiation for Buffalo and Dallas. The maximum and minimum ambient 
temperature and solar radiations occur with a time lag for these cities (Figure 5) where the interplay of energy exchange 
processes on the surface seems to favor the most conductive intermediate layer versus the least conductive layer or 
vice versa in certain months of the year. 

TDA was selected by the algorithm in hottest months of the year in Miami (table 8) to yield the highest energy 
dissipation rates. One interesting observation from the monthly results is that only the two intermediate material 
with lowest and highest thermal conductivities were selected as the choice blanket material through the optimization 
algorithm. The fact that TDA (the least conductive) or saturated clay (the most conductive) were selected among the 
provided list of intermediate material suggests that the ground pipes can benefit from a blanket above the pipe burial 
depth throughout the year but not necessarily from the one with highest insulation properties. An energy weighted 
average value of the operating parameters for heating and cooling seasons in selected cities is also provided at the end 
of each table for the sake of comparison between cities. 

It should be noted that this modeling procedure do not take into account other characteristics of either of the 
TDA or saturated clay materials. Characteristics such as the porosity and water holding capacity which might poten- 
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Table 7: GA results for Dallas 



Month 


Fluid 


EWT m i n 


EWT max 


EWT mean 


Thickness 


Position 


Intermediate 


Energy 


% 




(index) 


(°0 


(°C) 


(°C) 


Cm) 


Cm) 

\ Lll J 


material 


(Watt) 




uct 


i 
i 








u.z 


U.J 


1 1JA 




U.Z 


Nov 


Z 


-J.U 




J.Z 


n a 


1 7 
1 .z 


Saturated clay 




A Q 


1 Y , , . 

Dec 


c 

J 


-J.U 




^ 7 
O.Z 


u. / 


1 1 


Saturated clay 


1 OQA 


1 7 


Jan 


A 

4 


-J.U 




J.Z 




U. 1 


1 1JA 


1 O^/l 




rcu 


Z 


-J.U 




7 1 


U.J 


1 .J 


Saturated clay 


1017 


A Q 


Mar 


l 


-3.0 




5.2 


0.9 


0.8 


Saturated clay 


937 


6.8 


Apr 


l 


-2.1 




5.2 


0.4 


1.5 


Saturated clay 


825 


6.7 


May 


l 




53.6 


41.4 


0.2 


1.0 


TDA 


1500 


0.3 


Jun 


3 




53.4 


36.2 


0.2 


0.7 


TDA 


1500 


1.7 


Jul 


1 




51.8 


29.5 


0.5 


0.8 


TDA 


1500 


1.7 


Aug 


4 




51.6 


46.9 


0.5 


0.9 


TDA 


1500 


0.9 


Sep 


3 




56.6 


35.1 


0.2 


0.1 


TDA 


1500 


1.7 


Avg. Heating 




-2.9 




5.5 


0.7 


1 




1000 




Avg. Cooling 






54.3 


38.7 


0.3 


0.7 




1500 





tially contribute to considerably difFerent performance results than that only based on the heat conduction in the soil 
medium. TDA's porous structure can potentially enhance the moisture migration to the underlying layers of soil where 
higher moisture can contribute to higher thermal conductivities of the soil around the pipes. It is expected that this 
characteristics of TDA has a more substantial effect in the summer time specially in regions with less rainfall events. 
Authors are currently in the process of installation of a field experiment with TDA blanket and a control clayey section 
to verify other potential enhancements on GSHP performance by utilization of the TDA layer. 

A comparison between the working fluid minimum and average temperatures for the heating season in Buffalo 
and Dallas shows that relatively similar ranges are selected by the algorithm for optimal operation in cold season. On 
the other hand, as expected, the optimized maximum and minimum fluid temperatures in summer are comparatively 
higher for Miami and Dallas than corresponding values for Buffalo to ensure the highest energy dissipation to the 
ground. 

The values in the last column in tables 6 to 8 represent the percentage energy extraction/dissipation rate increase 
compared to the homogenous soil profile in different months. There are trends in comparing the results for similar 
seasons in different cities. The percentages increase in energy extraction rate for the coldest months in Buffalo (Jan- 
Feb) are as high as 18-19 % as compared to similar time periods for Dallas with highest values of 5-6 %. This finding 
can be translated in the potential for higher efficiency achievement with a non-homogeneous soil profile in heating 
season with more pronounced effects in colder climates (Buffalo). It should be noted again that the material choice 
for the intermediate layer in the cold season in Buffalo is TDA versus the saturated clay for Dallas. 

A similar comparison for the cooling season in all cities reveals an interesting observation. The highest efficiency 
improvements for cooling in Miami happen to be in the coldest month of the year with increase in energy dissipation 
rates of approximately 6-8 %. The optimization algorithm tends to choose higher working fluid temperatures, com- 
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Table 8: GA results for Miami 



Month 


Fluid 


EWT m i n 


EWT max 


EWT mean 


Thickness 


Position 


Intermediate 


Energy 


% 




(index) 


(°C) 


(°C) 


(°C) 


On) 

V AA V 


(m) 

V AA V 


material 


(Watt) 




uct 


z 




Jy.O 


44.0 


fl A 
U.O 


1 0 

1 .z 


Saturated Clay 


1 AQQ 


z.z 


Nov 


J 




1 Z.4 


J4.0 


U.3 


U.J 


Saturated Clay 


1JUU 


O 1 
U. 1 


1 Y , , . 

Dec 


Z 




7Q C 

/y.o 


AQ A 

oy.o 


U. / 


1 .u 


Saturated Clay 


1 

1 Jjy 


n 1 


Jan 


i 
1 




0 1 .J 


AT 7 
U / . J 


U.o 


1 1 
1 . 1 


Saturated Clay 


/ 1Z 


n a 
1 .4 


rcu 


i 
1 




oU.U 


AQ 9 

ov.z 


U.o 


1 .u 


Saturated Clay 


AzLQ 


A 1 
O.Z 


Mar 


1 




80.0 


69.2 


0.9 


1.0 


Saturated Clay 


1145 


8.5 


Apr 


4 




66.7 


66.0 


0.6 


1.0 


Saturated Clay 


1500 


0.6 


May 


1 




57.5 


44.6 


0.5 


0.1 


TDA 


1500 


1.1 


Jun 


2 




51.6 


50.3 


0.5 


0.7 


TDA 


1500 


0.2 


Jul 


3 




52.2 


51.3 


0.1 


0.3 


TDA 


1500 


0.9 


Aug 


3 




52.7 


43.1 


0.5 


0.1 


TDA 


1500 


0.6 


Sep 


1 




53.7 


51.7 


0.3 


0.1 


TDA 


1500 


0.3 


Avg. Cooling 






63 


55.2 


0.5 


0.6 




1322 





pared to the warmer months, with the highest intermediate layer conductivity (saturated clay) to maximize the heat 
flux to the ground and subsequently increase the efficiency. This trend does not repeat itself in warmer months as 
the algorithm searches through the large population of possible solutions. The comparison for the cooling season in 
all cities shows that relatively small efficiency improvements are achievable with the non-homogeneous soil profile 
for cooling purposes in warmer months of the year. It should be noted that the performance of the homogeneous 
system is very close to the maximum expected cooling capacity (1500 Watts) in warmest months of the year for all 
of these cities which contributes to the marginal improvements with a non-homogeneous soil profile. In other words, 
if the designed ground pipe network was supposed to provide higher energy rates, the percentage increase in energy 
extraction rates would probably be higher and more pronounced with the non-homogeneous soil profile. It would be 
safe to state that if chosen to use a non-homogeneous soil profile for cooling-only purposes in a warm climate, a high 
conductive intermediate layer can be more beneficial in coldest months of the year. 

After performing the analysis for each month and gaining adequate knowledge of the relationship between se- 
lected variables, the next step of the modeling was focused on finding the optimal values of the inlet fluid temperature 
and the configuration of the intermediate layer year-round. Knowing what the two dominant choices of the algorithm 
for the intermediate blanket are, separate runs of the model with TDA were executed to compare the results with the 
homogeneous and saturated clay soil profiles. The working fluid properties was set to 13 % Propylene-glycol and wa- 
ter mixture for all the annual runs. The existing knowledge of the system from monthly results helped constructing a 
more efficient optimization scheme with less variables for the annual runs. The goal for this step was to obtain the op- 
timal values of temperatures and intermediate layer configuration by optimizing the energy extraction/dissipation for 
the whole year. This approach provides the designer with the optimized working fluid temperatures and intermediate 
layer's configuration which yields the best performance. 

The results of annual runs are presented in tables 9 to 1 1 . The percentages in the table refer to the percent increase 
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in the energy extraction/dissipation rates compared to the homogeneous soil profile. The obtained optimized tem- 
perature values and intermediate layer configuration for TDA were used to run the homogeneous and saturated clay 
scenarios for the sake of comparison. For the Buffalo case, the algorithm yields higher maximum and average tem- 
perature values for heating season compared to average monthly values (Table 6) to obtain the maximum performance 
over the entire season. A similar trend was observed for Dallas with higher mean working fluid temperature in heating 
season compared to average monthly values (Table 7). The working fluid temperature values are comparatively lower 
for Buffalo than Dallas for annual results. The optimized annual TDA thickness and positions are relatively smaller 
than the values obtained in the monthly optimization. The average annual attainable energy rates are also consistently 
lower for the annual runs compared to monthly values. 



Table 9: Annual results for Buffalo 



FWT 

vv L mm 

CO 


E W T mean _ heating 

CO 


Thickness 
(m) 


Position 
(m) 


EWT max 

CO 


E W J mean-cooling 

CO 


-2.3 


7.0 


0.5 


0.3 


45.7 


31.7 


TDA energy 
(Watt) 


Clay energy 
(Watt) 


Homogeneous energy 
(Watt) 


Mode 


% TDA 


% Clay 


324 


274 


280 


Heating 


15.7 


-2.1 


1497 


1375 


1391 


Cooling 


7.6 


-1.2 


Table 10: Annual results for Dallas 


EWT m i n 
CO 


E W 7 mean - heating 

CO 


Thickness 
(m) 


Position 
(m) 


EWTmax 

CO 


EWT mean - coo ii n g 

CO 


-2.7 


23.5 


0.5 


0.4 


52.6 


29.6 


TDA energy 
(Watt) 


Clay energy 
(Watt) 


Homogeneous energy 
(Watt) 


Mode 


% TDA 


% Clay 


635 


608 


612 


Heating 


3.8 


-0.7 


1431 


1378 


1382 


Cooling 


3.5 


-0.3 


Table 1 1 : Annual results for Miami 


EWT m i n 
CO 


E W T mean - heat in g 

CO 


Thickness 
(m) 


Position 
(m) 


EWT max 

CO 


EWJ ''mean-cooling 

CO 






0.8 


0.4 


73.0 


55.0 


TDA energy 
(Watt) 


Clay energy 
(Watt) 


Homogeneous energy 
(Watt) 


Mode 


% TDA 


% Clay 


1298.0 


1482.0 


1456.0 


Cooling 


-10.9 


1.8 



The values of the annual energy rates with TDA and saturated clay compared to the homogeneous case suggest 
considerably lower annual performance with clay for both heating and cooling season in Buffalo and Dallas (Tables 
9 and 10). The results for Buffalo are in the order of 1-2 % worse than the homogeneous case with clay. The similar 
comparison for Dallas shows nearly identical annual performance with clay as the homogeneous case. TDA exhibits 
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higher performance in heating than cooling (15.7 versus 7.6), and comparatively higher increase in energy rates for 
Buffalo than Dallas in heating (15.6 versus 3.8). The results for Miami on the other hand suggest a significantly better 
performance with clay layer than TDA (1.8 % versus -10.9 %) for cooling. Although, the value of annual increase in 
energy dissipation rate (1.8 %) with clay is considerably lower than the monthly values obtained for Miami (Table 8). 

The comparison between the optimized annual heating/cooling energy rates with the monthly values shows that 
the optimization algorithm yields values that stray from the optimized monthly values. This is because of the fact that 
it is more challenging for the algorithm to find operational parameters which level off the difference between seasonal 
energy rates and highest achievable monthly values. This causes the annual energy rates to be relatively lower than the 
ones attainable if the monthly control on the operating parameters (working fluid temperatures and intermediate layer 
configuration) was an option. These findings suggest that utilization of control strategies which allow the heat pump's 
functionality with variable working fluid temperatures (or variable refrigerant flow schemes) in different seasons can 
improve the overall performance of the GSHP. 

6. Conclusions 

Performance of a GHSP system with a non-homogeneous soil profile was examined via GA optimization algo- 
rithm. The evolutionary algorithm was given a range of working fluid properties, intermediate layer thermal properties, 
a range of operating temperatures, and the intermediate layer configuration to search for the optimized condition of 
the system. The optimization process was performed for three cities representing a cold (Buffalo), moderate (Dallas) 
and warm (Mmiami) climate for the sake of comparison. A summary of findings are listed below: 

• A generic non-homogeneous soil profile showed the potential for significant GSHP performance improvement 

• TDA demonstrated higher benefits in colder climate and with higher magnitudes in heating season 

• TDA demonstrated marginal enhancement in cooling due to insignificant capacity difference than homogeneous 
scenario 

• Saturated clay demonstrated potential for efficiency improvement in cooling for colder months in a relatively 
warm climate 

• The annual optimized energy rates are lower than monthly values suggesting the necessity of a monthly opera- 
tion control strategy 

Despite different performance achievements with either a low conductance (TDA) or a high conductance (saturated 
clay) intermediate layer, a non-homogenous soil profile demonstrated the potential for increasing the the GSHPs 
performance. A shift in perspective toward control strategies for GSHPs from only the heat pump side to the ground 
pipe side of the system is suggested based on the model results. Further investigation of other attributes of a layered 
system which can potentially enhance the GSHP's performance is still required. TDA for example has a very porous 
structure enabling the water passage through the intermediate blanket. This characteristics of TDA has the potential to 
increase the thermal conductivity of the underlying soil by increasing the soil moisture content. Further investigation 
of the findings is planned through a field demonstration recently built on the University at Buffalo properties. 
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